Expt 3

library(janitor)

## 
## Attaching package: 'janitor'

## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(here)

## here() starts at /Users/jasongeller/Desktop/SF_Expt

library(afex)

## Loading required package: lme4

## Loading required package: Matrix

## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4

## ************
## Welcome to afex. For support visit: http://afex.singmann.science/

## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'KR', 'S', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - NEWS: library('emmeans') now needs to be called explicitly!
## - Get and set global package options with: afex_options()
## - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************

## 
## Attaching package: 'afex'

## The following object is masked from 'package:lme4':
## 
##     lmer

library(emmeans)
library(data.table)
library(patchwork)
library(data.table)
library(BayesFactor)

## Loading required package: coda

## ************
## Welcome to BayesFactor 0.9.12-4.2. If you have questions, please contact Richard Morey (richarddmorey@gmail.com).
## 
## Type BFManual() to open the manual.
## ************

library(tidyverse)

## ── Attaching packages ──────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.1
## ✓ tidyr   1.1.1     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0

## ── Conflicts ─────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::between()   masks data.table::between()
## x tidyr::expand()    masks Matrix::expand()
## x dplyr::filter()    masks stats::filter()
## x dplyr::first()     masks data.table::first()
## x dplyr::lag()       masks stats::lag()
## x dplyr::last()      masks data.table::last()
## x tidyr::pack()      masks Matrix::pack()
## x purrr::transpose() masks data.table::transpose()
## x tidyr::unpack()    masks Matrix::unpack()

library(ggplot2)
library(see)

setwd(here::here('Expt3_data', 'Gorilla_raw_data'))

data=here::here('Expt3_data', 'Gorilla_raw_data')  # path to data files

file_list=list.files(data, pattern=".csv") # list of data files
 
# read in all files
dataset <-
  do.call("rbind", lapply(file_list, FUN=function(files){
    
    for (i in 1:length(files)){ 
      if(file.exists(files[i])){
        message( "now processing:", files[i])
      }
    }
    fread(files, header=TRUE, sep=",", na.strings = "", fill=TRUE)})) #fread makes reading in files quick

## now processing:data_exp_15395-v4_task-7kze.csv

## now processing:data_exp_15395-v4_task-bqol.csv

## now processing:data_exp_15395-v4_task-u417.csv

## now processing:data_exp_15395-v4_task-wqdy.csv

## now processing:data_exp_15395-v5_task-7kze.csv

## now processing:data_exp_15395-v5_task-bqol.csv

## now processing:data_exp_15395-v5_task-u417.csv

## now processing:data_exp_15395-v5_task-wqdy.csv

#

dd<-dataset %>% 
  janitor::clean_names(.) %>%
  dplyr::filter(zone_type=="response_button_text")

#response as character
dd$response<-as.character(dd$response)

dd

#get Rts
rt<-dataset %>% janitor::clean_names(.) %>%  dplyr::filter(zone_type=="continue_button", display=="study") 

# get RT
rt$reaction_time<-as.numeric(rt$reaction_time)
rt1<- rt %>% 
  dplyr::group_by(participant_private_id , condition) %>% 
  dplyr::summarise(mean=mean(reaction_time, na.rm=TRUE))

## `summarise()` regrouping output by 'participant_private_id' (override with `.groups` argument)

rt2=rt1 %>% tidyr::pivot_wider(names_from=condition, values_from = "mean")

ttestBF(x=rt2$normal, y=rt2$SF, paired=TRUE, data=rt2)

## Warning: data coerced from tibble to data frame

## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 0.1569365 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS

rt2

#6.67

#response as character

ex3=dd %>% dplyr::mutate(condition1= dplyr::case_when( 
  condition == "SF" ~ "Sans Forgetica", 
  condition =="normal" ~  "Arial", 
), isold= dplyr::case_when (
  old_new== "old" ~ 1, 
  old_new== "new" ~ 0), 
sayold=dplyr::case_when( 
  response=="old"~ 1, 
  response=="new" ~ 0, 
  ))


#classic SDT for those wanting to compare
sdt <- ex3 %>% 
  dplyr::mutate(type = "hit",
         type = ifelse(isold==1 & sayold==0, "miss", type),
         type = ifelse(isold==0 & sayold==0, "cr", type),  # Correct rejection
         type = ifelse(isold==0 & sayold==1, "fa", type))  # False alarm

sdt <- sdt %>% 
  dplyr::group_by(participant_private_id, type, condition1) %>% 
  dplyr::summarise(count=n()) %>% 
  tidyr::spread(type, count)  # Format data to one row per person

## `summarise()` regrouping output by 'participant_private_id', 'type' (override with `.groups` argument)

sdt <- sdt %>% 
  dplyr::group_by(participant_private_id, condition1)%>%
  dplyr::mutate(hr = hit / (hit+miss),
         fa = fa / (fa+cr)) %>%
  dplyr::mutate(hr=dplyr::case_when(
    is.na(hr) ~ 0.99,
    TRUE ~ hr), 
    fa=dplyr::case_when(
      is.na(fa)  ~ 0.01,
    TRUE ~ fa),
     zhr=qnorm(hr), 
     zfa=qnorm(fa), 
    dprime = zhr-zfa) %>%
  ungroup()

sdt

# get prob of responding old 
sdt1=sdt  %>% dplyr::select(participant_private_id, condition1, hr, fa) %>% 
  tidyr::pivot_longer(hr:fa, names_to="type") %>%
  dplyr::mutate(isold=case_when(type=="hr" ~ "Old", type=="fa" ~ "New"))



sdt1$isold<-factor(sdt1$isold, levels=c("Old", "New"))

sdt1$Condition<-factor(sdt1$condition1, levels=c("Sans Forgetica", "Arial"))


dprimebf=sdt  %>% dplyr::select(participant_private_id, condition1, dprime) %>%
  tidyr::pivot_wider(names_from=condition1, values_from = dprime)

#paired ttest
a1 <- t.test(dprimebf$Arial,dprimebf$`Sans Forgetica`, paired=TRUE, data=dprimebf)

#Bayes Factor
dprimebf<-as.data.frame(dprimebf)

ttestBF(x=dprimebf$Arial, y=dprimebf$`Sans Forgetica`, paired=TRUE, data=dprimebf)

## Bayes factor analysis
## --------------
## [1] Alt., r=0.707 : 0.1466792 ±0%
## 
## Against denominator:
##   Null, mu = 0 
## ---
## Bayes factor type: BFoneSample, JZS

#BF=7.14

# path to data files

library(Rmisc)

## Loading required package: lattice

## Loading required package: plyr

## ------------------------------------------------------------------------------

## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)

## ------------------------------------------------------------------------------

## 
## Attaching package: 'plyr'

## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following object is masked from 'package:purrr':
## 
##     compact

## The following object is masked from 'package:here':
## 
##     here

means <- sdt1 %>% mutate(Condition=case_when(
  condition1=="SF"~ "Sans Forgetica", 
  condition1=="normal"~ "Arial"
))

means$Condition<-factor(means$condition1, levels=c("Sans Forgetica", "Arial"))



oldnewsub=summarySEwithin(data = means, measurevar = "value",
                       withinvars = c("Condition", "isold","type"), idvar = "participant_private_id")

## Automatically converting the following non-factors to factors: type

bold <- element_text(face = "bold", color = "black", size = 14) #axis bold

p1<- ggplot(means, aes(Condition, value, fill=Condition))+
  facet_grid(~isold) + 
  geom_violin() + 
  geom_jitter2(width=0.11, alpha=.5)+ 
  geom_line(data=oldnewsub,aes(y=value, group=1), size=1)+ 
  geom_pointrange(data=oldnewsub, aes(y=value, ymin=value-ci, ymax=value+ci), size=1, color="white")+ 
  theme_bw(base_size=14)+
  labs(y="Pr Saying Old", x="Font Type") + 
  theme(legend.position = "none") + 
  theme(axis.text=bold) 


#oldnew=brm(glmm2, data=ex3, family=bernoulli(link="identity"), prior=Priors, sample_prior = TRUE,  cores=6, inits = 0, control = list(adapt_delta = .9), iter=3000)

p1

library(ggplot2)

bold <- element_text(face = "bold", color = "black", size = 14) #axis bold

sdt$condition<-factor(sdt$condition1, levels=c("Sans Forgetica", "Arial"))


dsw=summarySEwithin(data = sdt, measurevar = "dprime",
                       withinvars = "condition1", idvar = "participant_private_id")

## Automatically converting the following non-factors to factors: condition1

p2<- ggplot(sdt, aes(condition1, dprime, fill=condition1))+
  geom_violin() + 
  geom_jitter2(width=0.11, alpha=.5)+ 
  geom_line(data=dsw,aes(y=dprime, group=1), size=1)+ 
  geom_pointrange(data=dsw, aes(y=dprime, ymin=dprime-ci, ymax=dprime+ci), size=1, color="white")+ 
  theme_bw(base_size=14)+
  labs(y="Sensitivity (d')", x="Font Type") + 
  theme(legend.position = "none") 

p2

patchwork1= p1/ p2 
patchwork1 + plot_annotation(tag_levels = 'A')

#Supplemental Analysis

#fit GLMM in brms to extract the BF

oldnewglme=glmer(sayold~isold*condition1+(1+condition1|Participant.Private.ID)+ (1+condition1|Stims), data=ex3, family=binomial(link="probit"))

prior<-prior(normal(0,1), class="b") # weakly informed prior on the coefficents

#fit the brms model
oldnewbrm=brm(sayold~isold*condition1+(1+isold*condition1|Participant.Private.ID)+ (1+isold*condition1|Stims), data=ex3, family=bernoulli(link="probit"), prior=prior, sample_prior = TRUE,  cores = 4) 


#get BF for interaction which is difference in dprime from brms model
dprime=hypothesis(oldnewbrm, 'isold:condition1 = 0')